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I. INTRODUCTION 



Understanding charge and energy transport at the nanoscale 
is essential for the design of stable and reproducible molecu- 
lar electronic components such as transistors, "refrigerators", 
and energy conversion devices IH]. While detailed model- 
ing is necessary for elucidating and optimizing the transport 
characteristics of such devices, in this paper we embrace an 
alternative-minimal approach [|2|]. With the motivation of ex- 
ploring the fundamentals of quantum transport in correlated 
electron systems, we focus on the dynamics of "impurity 
models" ||3[], consisting a small subsystem (molecule, quan- 
tum dot) interacting with two electronic reservoirs, driven to 
a nonequilibrium steady-state by a DC voltage bias. While 
the impurity object includes only few degrees of freedom, it 
incorporates many-body interactions, making exact analytical 
solutions generally inaccessible. Among the standard models 
considered in this context are the single impurity Anderson 
model (SIAM), combining a single electronic level with up 
to two interacting electrons coupled to metallic leads and 
the spinless two-level Anderson model (2LAM), consisting a 
spinless dot with two interacting (HOMO and LUMO) levels 
hybridized with electronic reservoirs [Hl-lzl. 

Even in the steady-state limit, the analysis of such nonequi- 
librium systems turns out to be intricate, and analytical so- 
lutions are lacking, see e.g., 1^. Various numerical simula- 
tion approaches have been developed, including perturbative 
treatments ^ and renormalization-group techniques llsl Uoll . 
Even more difficult is the description of the time evolution of 
the system from some initial preparation towards steady-state 
under a finite voltage-bias. The transient nonequilibrium dy- 
namics of the Anderson model, and its variants, has been re- 
centWsimulated using path-integral Monte-Carlo simulations 
lfTTl - fl3ll and influence-functional methods lITii ITsll . Several 
factors should be considered for fully understanding the dy- 
namics of such models: (i) the finite external bias, driving the 
system out-of-equilibrium, (ii) electron-electron interaction. 



or more generally many-body interactions, (iii) band- structure 
effects, and (iv) the device temperature. The combined effects 
of these four ingredients on the time evolution of a nanoscale 
object have not yet been fully understood ITgIi . 

Our objective here is to follow the dynamics of simple 
nanoscale junctions employing the SIAM and the 2LAM 
models as prototypes. We explore the role of the temper- 
ature and interaction strength in determining both the short 
time evolution of the system and its steady-state proper- 
ties. For achieving this task we adopt the recently developed 
numerically-exact influence functional path integral (INFPI) 
technique II15I1 . This method relies on the observation that in 
out-of-equilibrium (and finite temperature) cases bath correla- 
tions have a finite range, allowing for their truncation beyond a 
memory time dictated by the voltage-bias and the temperature. 
Taking advantage of this fact, an iterative-deterministic time- 
evolution scheme has been developed where convergence with 
respect to the memory length can in principle be reached. As 
convergence is facilitated at large bias, the method is well 
suited for the description of the real-time dynamics of single- 
molecule devices driven to a steady-state via interaction with 
biased leads. In this respect the INFPI approach is comple- 
mentary to methods applicable predominantly close to equi- 
librium, e.g., numerical renormalization group techniques yl. 

The principles of the INFPI approach have been detailed 
in Ref. il5il , where it has been adopted for investigating, at 
zero temperature, dissipation effects in the nonequilibrium 
spin-fermion model, and the population dynamics in a cor- 
related quantum dot, investigating the Anderson model. The 
focus of the present study are transport characteristics of cor- 
related nonequilibrium models, thus we introduce the INFPI 
approach in this context only. We demonstrate that the method 
can feasibility treat various impurity models. In particular, the 
population dynamics and the electron current in the SIAM and 
the 2LAM models are simulated at nonzero temperatures. 

The paper is organized as follows. In Sec. II we describe 
the INFPI method in the context of quantum transport junc- 
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tions. The nonequilibrium dynamics of the Anderson dot is 
studied in Sec. III. The spinless two-level Anderson model is 
discussed in Sec. IV. Some conclusions follow in Sec. V. 




FIG. 1: Schematic representation of the two models considered in 
this work: (top panel) the single impurity Anderson model with on- 
site repulsion tenns on the dot; (bottom panel) the spinless two-level 
Anderson model, with two electronic levels allowing for up to two 
interacting electrons. 



II. GENERAL FORMULATION 

We detail here the INFPI method in the context of quan- 
tum transport models. A more general presentation, dealing 
with both transport and dissipation in nonequilibrium open 
systems, is given in Ref. lITsil . The generic setup consid- 
ered includes a quantum impurity (subsystem) coupled to two 
metal leads (reservoirs) driven to a non-equilibrium steady- 
state through the application of a finite DC voltage-bias. The 
quantum impurity may be realized by a magnetic impurity, 
double quantum dots, or a multi-state quantum dot. The elec- 
trodes are modeled by two fermionic continua. System-bath 
couplings allow for a particle transfer between the impurity 
and the leads. We assume that the reservoirs' electrons are 
non-interacting, and include many-body interactions within 
the subsystem only, accounting for an additional energy cost 
for double occupancy. For a schematic representation see Fig. 
[U Our generic Hamiltonian is given by 

H = Ho + Hi, (1) 

where includes the exactly solvable noninteracting part 
combining the two leads, the noninteracting part of the sub- 
system, and impurity-bath hybridization terms. Many body 
interactions are incorporated into Hi, and we confine our 



present analysis to the special form 

Hi = U[nin2 - ^{ni + 712)]. (2) 

Here rii are occupation number operators for the subsystem 
with U as an interaction parameter. The states ' 1 ' and '2' may 
either symbolize the spin orientation, or count the (subsystem) 
electronic states. This structure allows for the elimination of 
Hi via the Hubbard-Stratonovich (HS) transformation UtI . 
In particular, in the Anderson model [see Eq. (fT7t l Hi ac- 
counts for the double occupancy energy cost on the dot. Sim- 
ilarly, in the 2LAM [Eq. ( |27] ll Hi constitutes the repulsion 
energy between electrons occupying the dot levels. 

Our objective here is to calculate the dynamics of a 
quadratic operator A, either given by subsystem or baths de- 
grees of freedom. This can be done by studying the Heisen- 
berg equation of motion of an exponential operator e^"^, with 
A a variable that is taken to vanish at the end of the calculation, 

(i(t)) = Tr(pi) = lim 4Tr[p(0)e'^*e^^e-^^*] . (3) 

A-s-O OA 

Here p is the total density matrix and the trace is performed 
over subsystem and reservoirs degrees of freedom. For sim- 
plicity, we assume that at the initial time {t = 0) the dot and 
the baths are decoupled, and that the baths are prepared in a 
nonequilibrium biased state. The time-zero total density ma- 
trix is therefore given by the product state p{0) = ps{0) <8> 
PL® PR- We proceed and factorize the time evolution operator 
using a standard breakup, e*^* ~ (e'^*')^, further assuming 
the Trotter decomposition e'™ w (^eiHoSt/2^iHiStpiHoSt/2y 
The many-body term Hi can be eliminated by introducing 
auxiliary Ising variables s = ± via the Hubbard-Stratonovich 
transformation iITtIi . 

±iHiSt _ }_\^ „-sK±{n2-ni) 



Here k± = k' =F in", k' = smh^^[siii{StU/2)]^/^, k" = 
sin-i[sin((5if7/2)]i/2. The uniqueness of this transformation 
requires USt < tt. In what follows we use the following short 
notation, 

g-H'±(s) = g-s/t-t(n2-ni) (-5-) 

Incorporating the Trotter decomposition and the HS transfor- 
mation into Eq. ([3]l, we find that the time evolution of A is 
dictated by 



p{0) ^e^f^oSt/2^iHiSt^iHoSt/2"^ gAA ^^-iHoSt/2 ^-iHiSt ^-iHoSt/2"^ 
lim A|J_ J dsfdsf...ds^Tr[p{0) (^e^HoSt/2^H+(s+,)^^HoSt/2^^ |^^^HoSt/2 ^H^(st) ^^HoSt/2^^ 



{sD -iHoSt/2 



g-i_f/o5t/2gff_(s-)g-t 



(6) 
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The above equation is exact in the limit St —>■ 0. We refer to where each element in the product, besides the first one, is 
the integrand as an "Influence Functional" (IF), and denote it given by the ratio between truncated IF, 
by /(s^, S2, sf ...s^). As discussed in Ref. 111511 . in standard 
nonequilibrium situations, even at zero temperature, bath cor- 
relations die exponentially, thus the IF can be truncated be- 
yond a memory time Tc = NgSt, corresponding to the time 

beyond which bath correlations may be controllably ignored. ^i^t ^ ■^fc+ii •*fc+Af,-i) 

Here Ng is an integer, and the correlation time Tc is dictated Is{sk, Sk+i, Sk+N^-i) — j^g± g± ^± y 
by the nonequilibrium situation, Tc ^ This argument ' ' '^+^=^2 

implies the following (non-unique) breakup lITsIl 

^-^s('S^_7V^_i_ij s^_7V,+2' '^w)' with 



-1) 



■Tr 



(9) 



Herea+(4) = [e'HoSt/2^H^(st)^HoSt/2\^ g ^ 
We now define the following multi-time object, 

and time-evolve it by multiplying it with the subsequent trun- 
cated IF, then summing over the intermediate variables, 

^(*fc+l' ■^fe+2' '^fc+iVs-l)-^.'5('^fe+l' •*fc+2' '^fc+jvj- 

(11) 

Summation over the internal variables results in the time local 
expectation value, e.g., at tk we get 

^ '-^) = ^ 7?.(Sj._,_2-Afsi '5fc+3_iV, ' )(12) 

This procedure is repeated for several values of small A. Tak- 
ing the numerical derivative with respect to A, the expectation 
value of the operator of interest, at a particular time, is re- 
trieved, {A{tk)). 

The truncated influence functional in Eq. (|9]l is the core 
of our calculation. Since it includes only quadratic operators 
ifisll . it can be exactly calculated utilizing the trace formula 
for fermions lfl9ll . 

Tr[e*^ie^^^..e^^''] = det[l + e"ie"^..e™'']. (13) 

Here irip is a single particle operator corresponding to a 
quadratic operator Mp = J2i ji^p)i,j'^i^j- 4 (Cj) 



fermionic creation (annihilation) operators. At zero temper- 
ature we can formally write Eq. (|9) as 

/ oc (0|e^^ie^^^..e^^''|0) = dct[e"ie"^..e™'']occ, (14) 

where |0) is the initial (zero temperature) state of the total sys- 
tem and the determinant is carried over occupied states only. 
At finite temperatures Eq. (|9]l can be represented by 

/ « TT[e^''e^'\..e^'-{pL ®Pr® PsiO))], (15) 

where pa corresponds to the time-zero density matrix of the 
a = L,R fermion bath. ps{0) denotes the subsystem initial 
density matrix. Assuming that these density operators can be 
written in an exponential form, e^^, with M a quadratic oper- 
ator llisll . application of the trace formula leads to 

= det{ [/l - h] «> [Ib - Ib] ® [Is - fs]. 

+ e™ie'"^..e"-[/L®/fl®/5]}. (16) 

The matrices and Is are the identity matrices for the a 
space and for the subsystem, respectively. The functions /l 
and ffj are the bands electrons' energy distribution, = 
jg;3„(£-MQ) _|„ 1]^^, with the chemical potential /i^ and tem- 
perature /3a- The subsystem (initial distribution) fs may vary, 
depending on the particular problem. For example, for the 
Anderson model (Sec. Ill) we consider a dot initially empty. 

In what follows we apply the INFPI method on two quan- 
tum impurity models, of interest in the context of molecular 
electronics, the SIAM and the 2LAM, see Fig. [T] with mini- 
mal modifications to the simulation code. Since both models 
admit the form ([Tli-(|2]i, we need only to separately construct 
the particular (single particle) noninteracting Hamiltonian Hq 
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and the operator of choice A. With this in hand, we can read- 
ily calculate the truncated IF of Eq. (|9]l and the ratio in Eq. (|8]l 
using the trace formula. We then time evolve the multi-time 
TZ structure following Eq. (flU . The time evolution of the 
operator of interest is acquired using Eq. (ITZi . We note that 
this iterative algorithm can be feasibly adopted for simulat- 
ing other models, including correlated multi-site chains with 
quartic interactions. However, the present implementation is 
limited by efficiency to models with two correlated sites ll20ll . 

Before discussing numerical results, we point out the dif- 
ferent sources of errors in our calculations, and explain how 
to control and overcome them. There are three sources of 
systematic error within our approach, (i) Bath discretiza- 
tion error. The electronic reservoirs are explicitly included 
in our simulations, and we use bands extending from — D to 
D with a finite number of states per bath per spin {Ls). This 
stands in contrast to standard approaches where a wide-band 
limit is assumed and analytical expressions for the reservoirs 
Green's functions are adopted iflll [l2l [Till . As we show be- 
low (see Fig. |6]l, by increasing the number of bath states Ls 
we can unequivocally reach convergence, typically employ- 
ing Ls > 100 states. We also note that while it is sometimes 
advantageous to encompass the leads' effect into self ener- 
gies terms, complex dispersion relations can be easily handled 
within our method, (ii) Trotter error. The time discretization 
error, order of {U5t)'^, originates from the approximate fac- 
torization of the total Hamiltonian into the non-commuting 
Hq (two-body) and Hi (many-body) terms, see text after Eq. 
(|3]l. While for J7 — >■ and for small time-steps 5t ^ Q 
the decomposition is exactly satisfied, for large U one should 
go to a sufficiently small time-step in order to avoid signif- 
icant error buildup. Extrapolation to the limit (5t — > is 
straightforward in principle lITsll . (iii) Memory error. Our 
approach assumes that bath correlations exponentially decay 
resulting from the nonequilibrium condition A/i ^ 0. Based 
on this crucial element, the influence functional may be trun- 
cated to include only a finite number of fictitious spins Ns, 
where Tc = Ns5t ^ 1/A/i for the population dynamics and 
Tc =^ 2/ A/^ for the particle current (see Figs. Island [TTl. The 
total IF is retrieved by taking the limit Ns N,(N ^ t/6t). 
However, one should be careful at this point: Increasing the 
memory length Tc by adding more and more Trotter-terms into 
the truncated IF [Eq. (|9]l] results in a build-up of the time 
discretization error, unless the time-step is controlled concur- 
rently. Thus, one should carefully monitor both the time-step 
and the memory size for achieving reliable results. This chal- 
lenge is similar to that encountered in the standard QUAPI 
method 1I21II22I1 . 

It should be noted that the convergence with respect to 
memory error is currently the most challenging aspect of the 
calculations with the INFPI approach. This limits us to rel- 
atively small values of the ratio of on-site correlation to hy- 
bridization strength. Future work will be devoted to algorith- 
mic optimization of the approach so that significantly larger 
memory times may be reached. 



III. ANDERSON DOT 

A. Model and Observables 

The single impurity Anderson Model (SIAM) [Ii]] is one 
of the most important models in condensed matter physics. 
While it was originally introduced to describe the behavior of 
magnetic impurities in non-magnetic hosts 1I23I1 . it has more 
recently served as a generic model for understanding quantum 
transport in correlated nanoscale systems ||24| - |26| . In such 
cases, the impurity is hybridized with two reservoirs main- 
tained at different chemical potentials, leading to nonequilib- 
rium particle transport. The model includes a resonant level 
of energy e^, described by the creation operator dj. (a =t, i 
denotes the spin orientation) coupled to two fermionic leads 
{a = L, R) of different chemical potentials /Iq,, but equal 
temperatures The Hamiltonian H — Hq + Hi [see Eqs. 
([T]l-(|2|l] includes the following terms 

+ ^^^'^^i,k.ada + h.C. 

Hi = U[nd,^nd,i ~ ■^{ud^^ + Ud^i)]. (17) 

Here ^ ^ ica,k,a) denotes the creation (annihilation) of an 
electron with momentum k and spin a in the a lead, U stands 
for the onsite repulsion energy, and Va^k are the impurity-a 
lead coupling elements. Ud.a = d-l-dcr is the impurity occu- 
pation number operator The shifted single-particle energies 
are denoted hy Ed = ed + U /2. We also define F — 
where Ta = ""X^fe |Kt,feP<5(e — £*.•) is the hybridization en- 
ergy of the resonant level with the a metal. In what follows 
we focus on two observables: the time dependent occupation 
of the resonant level and the tunneling current through the dot. 
The population dynamics {iid.a (t)) can be obtained by substi- 
tuting 

A = nd^a (18) 

in Eq. (|6]l. The current at the a contact (la.a) may be resolved 
in two ways. We may either calculate the population depletion 
(or gain) in the a lead by defining A as the sum over the a- 
bath number operators, 

^ = E<'o..Co,.,.. (19) 

k 

The current itself is given by the time derivative of the A ex- 
pectation value, {Ia,a) ~ ^(^(^))- Alternatively, the current 
at each end can be directly gathered by adopting the expres- 
sion A = — 25^j, VaMcl^ k a'^'^' "^i'h 3 as the imaginary 
part. In practice, we have employed the symmetric definition 

A= -^Y.^LA,k,ad. +^Y.^B^-^''kk,ad.. (20) 
k k 
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since its expectation value directly produces the symmetrized 
current 



(21) 




FIG. 2: Population of the resonant level in the Anderson model U = 
(thick full), U = 0.1 (dashed), U = 0.3 (dashed-dotted), U = 0.5 
(dotted). The physical parameters of the model are _D = 1, A/i = 
0.4, Ed = 0.3, rc«=0.025, and fiF = 10. The numerical parameters 
used are Ls = 240 lead states, Tc = 3.2 with Ns = 4 and 5t = 0.8. 
The U = case is compared to the wide flat band limit, Eq. i22\ 
(thin full line). 



B. Results 

We focus on the following set of parameters: a symmetri- 
cally distributed voltage bias between two leads with A/i ~ 
0.4, flat bands centered at zero (the Fermi energy) with a cut- 
off at _D = ±1, a resonant level energy Ed = 0.3, a hybridiza- 
tion strength Fq, = 0.025 = TT\Va.k\'^ Pa, with a constant den- 
sity of states Pa, onsite repulsion U/T ^2 — 10, and a zero 
magnetic field. For these parameters a convergence analysis 
carried out in Ref. fl3\ has revealed that supplying Lg > 100 
states per spin per bath suffices for mimicking a continuous 
band structure. We have also found that for Ap, = 0.4 a mem- 
ory size Tc '~ 1/A/i ^ 3.2 has lead to the convergence of the 
dot occupation when St = 0.8 and Ns — 4, provided ^ < 3 
ifTsi I27I1 . As we show below, the simulation of the current 
turns out to be more challenging as a larger memory size is 
required for reaching converging behavior, Tc ^ 2/Ap. 

Before presenting our results we clarify the initial condi- 
tions adopted here. As explained above, at i = we set the 
reservoirs and the system in a factorized state: The dot is as- 
sumed to be empty, and the two reservoirs are decoupled, each 
maintained in a canonical state characterized by the Fermi- 
Dirac statistics. This scenario is distinct from the interaction 
and voltage quenches considered in Ref. ifot] . 

Fig. |2] displays the time evolution of the dot occupancy 
{nd.cr) with increasing on-site interaction for fiV = 10, es- 
sentially reproducing the T = data of Ref. ifTsll . Details 
about convergence issues, and a comparison to Monte-Carlo 
data were included in Ref. ITsIIztIi . In order to examine the 
effect of the bandwidth on the details of the dynamics the evo- 
lution of the noninteracting case {U = 0) is further compared 
to the wide flat band (WFB) behavior lfl6l|28ll . 
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FIG. 3: Population of the resonant level in the Anderson model [/ = 
(thick full), U = 0.1 (dashed), U = 0.3 (dashed-dotted) at various 
temperatures, PV = 1, 2.5 and 10; top to bottom. Other parameters 
are the same as in Fig. |2] The inset compares the long-time U = Q.l 
behavior (□) to mean-field results (o) obtained from Eq. ( 123) . 



{nd.cr{t)) 



de[fL{e) + fR{e)] 



27r 

1 + e^^rt _ 2e-^* cos[(e - ed)t] 
r2 + (e - edY ■ 



(22) 



We find that the D/T = 2Q case inspected here deviates from 
the WFB result in both the short time behavior and the long 
time characteristics. However, general trends are maintained. 
We have also verified (data not shown) that the INFPI results 
approach the WFB limit when increasing the bandwidth, for 
U = 0. 

The effect of the temperature at different interaction 
strengths is analyzed in Fig. |3] adopting fiT = 0.1 — 10. For 
L/ = 0.1, a comparison between the long time INFPI limit and 
the mean-field theory HH, 



(rid,o-(i 00)) 



F 
2^ 



/L(e)+/fl(e) 



de. 



(23) 



reveals a good agreement (inset, Fig|3]l. 

For the same set of parameters we calculate next the sym- 
metric tunneling current {Ia{t)) through the SIAM. Simula- 
tion results for U = and U = 0.1 are presented in Fig. |4] 
The current enhancement with U can be reasoned by noting 
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FIG. 4: Current through the Anderson dot, U = (small dots), U = 
0.1 (large dots), Ed = 0.3, T = 0.05, /ST = 10. The [/ = case 
is compared to the WFB limit obtained from Eq. ll24t (thin full line). 
The numeric parameters are St — 1.6, Ns — 5 and Ls — 120. 
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The convergence of the tunneling current with respect to the 
number of bath states, time-step, and memory size has been 
carefully tested. In particular. Fig. ^demonstrates the behav- 
ior of the current with increasing memory size Tc = NgSt, 
showing that convergence is reached when Tc ~ 7 — 8. We 
note that a significantly shorter memory size (tc ^ 3 — 4) has 
been required for converging the dot occupancy ijTsj. This 
difference could be reasoned as follows. Since the tunnel- 
ing current is calculated at a specific contact, the memory size 
that should be accounted for inside the influence functional ^ 
should roughly scales with the bias difference at that contact. 
Thus, ^ A/i/2. In contrast, the population dynamics 
is sensitive to the full bias drop A/i, therefore bath correla- 
tions can be safely truncated beyond ^ 1/A/i. In Fig. 
|6] we present the behavior of the current upon increasing the 
number of bath states. It is interesting to note that the choice 
Lg = 40 states per spin per bath already reproduces results 
in a good agreement with the — s- oo limit. Thus, the finite 
temperature algorithm adopted here [Eq. (|T6^1. is superior 
to the strictly zero temperature algorithm of Ref. ifisll . even 
when applied to relatively low temperatures. 

It is also of interest to examine the temperature dependence 
of the asymptotic electric current. This information is con- 
veyed in Fig. |2]for zero and finite U using data at tT ~ 5. 
Results are also compared to the mean-field wide-band ap- 
proximation ll29i[30ll . 
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rVL{e)-fR{e)] 
{e-ed-U{nd,-a)f +r- 



;de. (26) 



FIG. 5; Convergence of the current {laii)) through the Anderson 
dot with increasing memory size Tc — NsSt. Ed = 0.3, U = 0.1, 
r = 0.05, /3r = 10. The numerical parameters are Ls = 120 states 
and St = 1.6. Ns ^ 2 (o), iV, = 3 (o), iV, = 4 (+), iV^ = 5 
(x), Ns — 6 (□), Ns = 7 (dotted line). Inset: zooming over the 
long-time values. 



Deviations from this result, for U = 0, indicate on the de- 
parture from the WFB approximation. In the large bias limit 
examined here (A/x/F — 8) the current saturates at low tem- 
peratures, /?r < 2.5, in agreement with the results of Ref. 



that the parameter Ed = ed + U/2 is fixed, thus the actual dot 
energy is down-shifted when increasing the interaction U. We 
again compare the noninteracting behavior with the dynamics 
in the WFB limit HI], 
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re-^*[/L(e 
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'(e-erf)2 + r2 
X {re~^VL{e) + fR{e)] 
^ Tcos[{e-ed)t][2fRie) + l] 
- (e-ed)sin[(e-erf)t][2/L(e) 

with the asymptotic value 

r hie) -Me) 



1]}, (24) 



{lL,a{t -> OO)) = 



27r 



r2 



(25) 



and {I]i^a-{t)) = — (/L.cr(— A/i, t)). Good agreement is ob- 
served in the long time limit. 
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FIG. 6; Convergence of the current (/^(t)) through the Anderson 
dot with increasing number of bath states Ls- Ed=0.3, U = 0.1, 
/3r = 10, F = 0.05, St = 1.6, iV. = 5. Ls=40 (heavy full), 80 
(dashed), 120 (dotted), 160 (dashed-dotted), and 240 (light full). The 
data lines for Ls > 80 are almost overlapping, see also the bottom 
inset. Top inset: Data as a function of Ls at Tt = 2.5. 
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FIG. 7: Steady-state current (/ss.o-) = {icr{t — >■ oo)) through the 
Anderson dot, (7 = 0.1 (o) and [/ = (□), Ed = 0.3, F = 0.05. 
The full lines are the results of a mean-field calculation, Eq. 
The numerical parameters are Ls = 120, Ns — 5 and 6t = 1.6. 



IV. SPINLESS TWO-LEVEL ANDERSON MODEL 

A. Model and observables 

The spinless two-level Anderson model (2LAM) and its 
extensions have been extensively studied in the context of 
molecular electronics, for exploring various effects in molec- 
ular conduction: vibrational effects llsill . thermoelectricity 
in molecular junctions lf32l fssll . radiation field-induced pro- 
cesses 1I34I1 . and Coulomb interaction effects Js]] . More re- 
cently, the mechanism of population inversion |5|] has been ex- 
plored using the asymmetric interacting 2LAM, where the two 
levels differently couple to the leads. Furthermore, by includ- 
ing a left-right asymmetry in the dot-leads coupling, the mech- 
anism of the transmission phase lapses in quantum dots llssll 
has been resolved within mean-field theories ll36[|37ll . Monte- 
Carlo techniques ITJ, and functional and numerical renormal- 
ization group approaches llssi l39ll . The 2LAM model incor- 
porates an impurity with two electronic levels ei < 62, de- 
scribed by the creation operator c?J,j, (m = 1, 2), coupled to 
two metal leads {a = L, R) of different chemical potentials. 
The Hamiltonian H = Hq + Hi includes the following terms 

Ho = (ei + C//2)ni + (£2 + C//2)7i2 + ^ 



a,A;,?7i— 1,2 



Hi 



Here 



U[nin2 - -{ni + 712)]. 



(27) 
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denotes the creation (annihilation) of an electron 

with momentum k in the a lead, ~ dl-^dm is the number 
operator for the impurity levels, and U is the charging energy. 
We also define the hybridization strength r„i = „i + F i}.,„ 
with Ta.m = TT J^k \^a.k.m\'^S{e — efc) and use flat bands ex- 
tending symmetrically between ±D. The dot shifted energies 
are denoted by i?,„ = Cm + U /2. This model is closely related 
to the interacting Anderson model analyzed in Sec. IH, tak- 
ing the two states here to emulate different spin orientations. 



However, here (i) only a single spin specie is considered, al- 
lowing for interference effects between the two transmission 
pathways, (ii) the dot levels are nondegenerate, and (iii) the 
impurity states differently couple to the leads, typically as- 
suming that the HOMO level, a deep molecular orbital, is cou- 
pled more weakly to the leads. 

The population dynamics of each electronic level (n,„(t)) 
and the current through the 2LAM are calculated numerically 
using the INFPI method, as prescribed in Sec. IF The current 
plotted will be the total symmetrized current flowing through 
the system, obtained by defining the operator of interest A as 



A 



k.m 



1 

0.8 
0.6 
0.4 
0.2 

0. 



^'^R,k,m"'^ 



(28) 
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FIG. 8: Population of the 2LAM electronic levels with increasing U 
term. (7 = (full), U = 0.1 (dashed), U = 0.2 (dashed-dotted). 
El = -0.1, E2 = 0.3, Fi.c = 0.025, T2,o. = 0.05, [3 = 200. The 
numerical parameters are 5t = 0.8, N3 = 5 and Ls = 120. 
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FIG. 9: Steady-state population of the 2LAM electronic levels: 
Comparison between the INFPI asymptotic data, extracted from Fig. 
[8](o), and mean-field results (□). 



B. Results 

We focus on the symmetric (L — R) case, and use 
the following set of parameters: FL.i=Fji' i=0.025 and 
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-U=0, (5=200 
■-U=0.1, p=200 
-U=0, P=20 
-U=0.1, P=20 

2 2^5 



FIG. 10: Current dynamics in the 2LAM with increasing U term. 
U = (full), (/ = 0.1 (dashed) and l3 = 200 (heavy) = 20 
(light). El = -0.1, E2 = 0.3, ri,„ = 0.025, Ta.c = 0.05. The 
numerical parameters are St — 0.8, Ns = 7 and Ls = 120. 
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FIG. 11: Convergence of the steady state current with increasing 
memory size Tc using different time steps St — 1.6 (empty sym- 
bols) and 5t=0.8 (full symbols) for P = 200 (circle) and /? = 20 
(square). Parameters are the same as in Fig. [To] 



the level's population satisfy 

{rimit 00)) 



r 

2^ 



U{nrn)r+Tl 



(29) 



where m = 2, 1 if m = 1, 2. In Fig. |9]we plot the asymptotic 
population dynamics, using the data from Fig. |8] and com- 
pare those values to mean-field results. As expected, the dis- 
crepancy between these two calculations increases for larger 
U. Deviations at [/ = probably stem from the fact that the 
INFPI method assumes finite bands of D ^ ±1, while mean- 
field results are calculated for WFB leads. 

We examine the temporal behavior of the current in Fig. 
[TOl varying the temperature and the many-body interaction 
strength. For the present set of parameters we conclude that 
the current decreases for large U, and that the temporal os- 
cillations are washed out with increasing temperature. Fi- 
nally, we use this data as highlighted in Fig. [TT]to expose 
a subtle convergence issue: the counteracting effect of differ- 
ent sources of errors, the time-step and the memoiy-size, and 
the challenge to overcome them both together. Employing the 
same set of parameters as in Fig. [TO] we extract the steady- 
state value for the current, and display it as a function of Tc, 
at two different temperatures, using two different time-steps. 
We find that for 4.5 < Tc < 8 the steady state results are 
almost fixed, fluctuating by only 1%. However, for Tc > 8 
a departure from the apparent steady state occurs, becoming 
larger for larger Tc- This behavior is caused by buildup of the 
Trotter factorization error within the truncated IF, Eq. (0. As 
expected, the error increases at larger U. To control this error, 
at large Tc a shorted time-step should be selected. 

Future work will be dedicated to the strong coupling limit, 
Tm > £2 ~ ci, for analyzing the charge oscillation effect IH]. 
The asymmetric L — R setup is also of great importance, for 
studying the phase lapses mechanism beyond the mean-field 
approximation, at strong driving ll36l . 



V. SUMMARY 



rL,2=rH,2=0.05. The bias (A^ = 0.4) will be symmetri- 
cally distributed between the leads, assuming flat bands cen- 
tered around zero with a cutoff at D = ±1. The interaction 
strength will be limited to U /Fi < 4 and the temperature 
will be varied between /3Ti ^ 1 — 10. Fig. [8] displays the 
levels' occupation as a function of time, for several interac- 
tion values, U=0, 0.1 and 0.2. We find that the HOMO pop- 
ulation {ni{t)) is increasing with U. In conjunction, due to 
the increased importance of repulsion effects on the dot, the 
LUMO population (n2(t)) depletes with U. The convergence 
of the data with respect to the number of bath states Ls, time 
step St, and memory size Tc ~ NgSt has been verified. A 
comparison to the WFB limit for the noninteracting case re- 
veals dynamical properties similar to those identified in Fig. 
121 Steady-state mean-field results are obtained by using ex- 
pressions analogous to Eqs. (|23] | and ( |26] | isj]. For example. 



We have employed here the INFPI method lITsIl for study- 
ing the population dynamics and the current behavior of two 
eminent molecular junction models: the single impurity An- 
derson model, and the 2-level Anderson dot. Considering 
voltage-biased junctions, the effect of the intra-dot electron- 
electron repulsion energy and the temperature were jointly 
analyzed. We have compared our results to mean-field calcu- 
lations, showing an increased discrepancy when many-body 
interactions are enhanced. A careful convergence analysis has 
been performed, demonstrating how to adequately converge 
the INFPI simulations. 

The INFPI method has been described here in connection 
with molecular transport junctions. We expect this flexi- 
ble tool to become useful for studying other-related impurity 
models, and for exploring nonlinear thermoelectric effects in 
molecular junctions llssll . In particular, future work will be fo- 
cused on simulating the dynamics of extended junctions, e.g.. 
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a multi-site chain, and on extending the method to include vi- like to acknowledge the NSF for financial support, 
brational effects Olll in a non perturbative manner. 
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